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Abstract 

We present a linear scaling formulation for the solution of the all-electron Coulomb 
problem in crystalline solids. The resulting method is systematically improvable and 
well suited to large-scale quantum mechanical calculations in which the Coulomb po- 
tential and energy of a continuous electronic density and singular nuclear density are 
required. Linear scaling is achieved by introducing smooth, strictly local neutralizing 
densities to render nuclear interactions strictly local, and solving the remaining neu- 
tral Poisson problem for the electrons in real space. While the formulation includes 
singular nuclear potentials without smearing approximations, the required Poisson so- 
lution is in Sobolev space H , as required for convergence in the energy norm. We 
employ enriched finite elements, with enrichments from isolated atom solutions, for an 
efficient solution of the resulting Poisson problem in the interacting solid. We demon- 
strate the accuracy and convergence of the approach by direct comparison to standard 
Ewald sums for a lattice of point charges, and demonstrate the accuracy in all-electron 
quantum mechanical calculations with an application to crystalline diamond. 

1 Introduction 

The evaluation of the electrostatic potential and total energy of crystalline solids has been an 
ongoing problem since the earliest days of solid state physics (|Madelung 1918[ lEwald 1921( 



Wiener and Seitz 1933j IFuchs 19351 llhm et al. 19791 IWeinert 198T]) . In ab initio density- 



functional ( Hohenberg and Kohn 1964 IKohn and Sham 19651 IJones and Gunnarsson 19891) 



calculations, the electrostatic (Coulomb) potential is a sum of nuclear and electronic con- 
tributions. In an infinite crystal, however, each of these terms diverges and the sum is 
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only conditionally convergent due to the long-range 1/r nature of the Coulomb interac- 
tion. Similarly, the total Coulomb energy is a sum of electron-nucleus, electron-electron, 
and nucleus-nucleus contributions, each of which diverges in an infinite crystal but combine 
to yield a finite total electrostatic energy per unit cell. Hence, in the all-electron quantum 
mechanical problem in solids, there are three distinct divergences which must be addressed 
simultaneously: (1) the 1/r divergence of the electrostatic potential at the nuclei, (2) the 
divergence of both potential and energy lattice sums due to the the long-range 1/r nature 
of the Coulomb interaction, and (3) the infinite self energies of the nuclei. 

It has been appreciated for some time that the divergences and conditional conver- 
gence of such extended lattice summations can be eliminated by formulating the summa- 
tions in terms of neutral densities that are well localized in real and/or reciprocal (Fourier) 
space (lEwald 19211) . In the pseudopotential approximation (IPickett 19891) . the nuclei and 
core electrons are replaced by smooth ionic cores, thus eliminating 1/r divergences and infi- 
nite self energies. In the conventional reciprocal space approach for such calculations in crys- 
tals (llhm et al. 1979( IPickett 1989|) . the remaining lattice sum divergences are eliminated by 
adding neutralizing densities to otherwise divergent Coulomb terms in such a way that the ef- 
fects of the added densities cancel in the final expressions. Remaining long-range interactions 
are then rendered short-ranged by transforming to reciprocal space, where smooth periodic 
functions, of infinite extent in real space, are well localized. However, the resulting expres- 
sions for the electrostatic potential and total energy contain structure factors and/or Ewald 
sums, and require at least 0(N log N) operations to evaluate, where N is the number of atoms 
in the unit cell. Moreover, since the approach relies on Fourier transforms, it is difficult to 
implement efficiently on large-scale parallel computational architectures due to the need for 
extensive interprocessor communications. The limitations of the reciprocal space approach 
have inspired much research on real-space and local-orbital based approaches flArias 1999[ 
IBeck 20001 ISoler et al. 20021 IPask and Sterne 2005at ITorsti et al."2006jl . which allow for 
better scaling, a variety of boundary conditions, and eliminate the need for Fourier trans- 
forms. These approaches accomplish a linear scaling solution of the Coulomb problem by 
solving Poisson's equation in real space, or evaluating the associated integral, thus alleviating 
the need for Fourier transforms, and allowing the use of efficient multi-level precondition- 
ing (IBrandt 19771 IBeck 20001) . 

In the all-electron quantum mechanical context, however, the divergences of the nu- 
clear potentials and self-energies must be confronted in addition to the divergences in lat- 
tice sums. Moreover, the rapid, local variation of the electronic density and potential 
in the vicinity of the nuclear singularities must be addressed. One approach to dealing 
with nuclear singularities is to approximate the singular nuclear densities by finite, local- 
ized functions, e.g., Gaussians or step functions (IMerrick et al. 19951 IModine et al. 19971 
IGoedecker and Ivanov 1998| Wang and Beck 2000 Suryanarayana et al. 2010). This makes 
possible the solution of a nonsingular total electrostatic potential, due to both nuclei and 
electrons, in a single linearly scaling step via solution of Poisson's equation with nonsingu- 
lar total (nuclear + electronic) density as source term. Furthermore, it makes possible the 
direct evaluation of the total Coulomb energy, with finite nuclear self-energy that is readily 
removed. However, as the widths of the model nuclear densities are decreased toward physi- 
cal nuclear dimensions (on the order of 10~ 5 a.u.) to achieve convergence, the potential in the 
vicinity of the nuclei and nuclear self-energies become correspondingly large, causing greater 



2 



absolute errors, which hinder the computation of accurate energy differences for different 
configurations. Moreover, finer resolution, and correspondingly more degrees of freedom, 
are required to resolve the more rapidly varying potential in the vicinity of the nuclei, thus 
increasing computational cost. This was well demonstrated in a remarkable calculation of 
the Coulomb potential of a uranium dimer (IGoedecker and Ivanov 1998|) using a second gen- 
eration interpolating wavelet basis, wherein 22 levels of refinement were required to reduce 
the maximum error in the potential to order 10 -2 a.u. These difficulties are a consequence 
of the divergent limit of the sequence: since the three-dimensional Dirac-delta is not in H^ 1 , 
the corresponding solution of Poisson's equation is not in H 1 , where H m is the Sobolev space 
of order m. Therefore, the Coulomb energy diverges and a convergent approximation in H 1 
cannot be constructed. 

Other approaches for all-electron quantum mechanical calculations include singular nu- 
clear contributions analytically and compute the remaining contributions analytically or 
numerically, depending on the choice of basis. One such approach is to compute the elec- 
tronic contributions via solution of Poisson's equation with nonsingular, though rapidly vary- 
ing, electronic density as source term (IWhite et al. 1989] IMurakami et al. 1992| Tsuchida 
and Tsukada IBatcho 20001 IBeck 20001 |Yamakawa and Hyodo 2005] ITorsti et al. ~2006l 



Bylaska et al. 2009 ILehtovaara et al. 2009|) ; where in crystalline calculations (Tsuchida and 



Tsukada I1995( IBeck 200~0~] ITorsti et al. 2006]) , the density must be neutralized to avoid di- 
vergent lattice sums. With mesh refinement (IMurakami et al. 19921 Tsuchida and Tsukada 
IT9951 IBatcho 20001 ITorsti et al. ~2006l |Bylaska et al. 2009| ILehtovaara et al. 2009j) and/or 
addition of well-chosen localized functions to the basis (Yamakawa and Hyodo 2005), one can 
then solve the resulting Poisson equation for the electronic contribution to the all-electron 
potential in a single O(N) scaling step. Integrals involving the singular total (electronic + 
nuclear) potential can be efficiently computed using a transformation of the singular nuclear 
part (IWhite et al. 19891 IMurakami et al. 19921 IBatcho 20001 |Yamakawa and Hyodo 2"005| ). 
However, in the crystalline case (ITsuchida and Tsukada 1995[) . the remaining nuclear contri- 
bution must then be computed by lattice summation, which scales as 0(N 2 ) or 0(iVlog N) at 
best. An alternative approach to the calculation of the electronic contribution is direct eval- 
uation of the associated Coulomb law integral (IGenovese et al. 2006| Juselius and Sundholm 
I2007| IWatson and Hirao 2008| ILosilla et al. 2010]) . This approach, recently extended to pe- 
riodic calculations ( ILosilla et al. 2010]) . can accommodate a variety of boundary conditions, 
can attain high accuracy, and when combined with the fast multipole method (Greengard 
and Rokhlin I1987| IStrain et al. 1996]) for far-field contributions, can achieve O(N) scal- 
ing (IWatson and Hirao 2008|) . However, the approach can be sensitive to the approximation 
employed for the singular l/|x — x'| kernel ( IJuselius and Sundholm 2007| Watson and Hi- 
rao |2008] ILosilla et al. 2010]) . and as with all such approaches computing just the electronic 
contribution, it leaves the singular nuclear contribution to be computed by other means. 

Another approach, employed in accurate all-electron density-functional electronic struc- 
ture methods ( Singh and Nordstrom 2006 ), employs a dual representation of the density and 
potential ( |Rudge 1969] IWeinert 198T1 IBlochl 19941 |Nikolaev and Dyachkov 12002] ) . In this 
approach, the unit cell is partitioned into atom-centered sphere and interstitial (between 
spheres) regions. Inside the spheres, where the potential is singular and most rapidly vary- 
ing, a radial- angular spherical harmonic representation is employed. Outside the spheres, 
where the potential is generally smooth, a Fourier representation is used. The all-electron 
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potential is then computed in two steps: (1) compute the smooth interstitial potential by 
standard Fourier techniques, and (2) solve the boundary value problem in each sphere using 
Green's functions, with boundary values from the previously computed interstitial poten- 
tial. By virtue of the dual representation, this method can provide highly accurate and 
efficient solutions. However, due to the reliance on Fourier expansions, the scaling is at best 
O(NlogN) and efficient parallel implementation is difficult due to the need for extensive 
interprocessor communications. 

Fast multipole methods (Greengard and Rokhlin 1987 IStrain et al. 1996 ) provide an 
O(N) scaling solution of the Coulomb problem and can accommodate periodic (Challacombe 
et al. [T9971 IKudin and Scuseria 19981 IKudin and Scuseria 20001 IKudin and Scuseria 2004")) 
as well as non-periodic boundary conditions. In the context of a Gaussian representation 
of the charge density, these have become the method of choice for large-scale calculations, 
as they can provide high accuracy, with linear scaling and a moderate prefactor, and are 
well suited to parallel implementation. However, the computational cost increases rapidly 
for higher-quality basis sets ( |Kurashige et al. 2007[ ), and due to the reliance on a Gaussian 
representation, it is not a general purpose method fIGenovese et al. 2006)1 . 

In this paper, we present a systematically improvable, linear scaling formulation for the 
solution of the all-electron Coulomb problem in crystalline solids. Linear scaling is achieved 
by introducing smooth, strictly local neutralizing densities to render nuclear interactions 
strictly local, and solving the remaining neutral Poisson problem for the electrons in real 
space. In so doing, the all-electron problem is decomposed into analytic strictly-local nu- 
clear, and numerical long-range electronic parts; with required numerical solution in H 1 
so that convergence is assured and approximation is optimal in the relevant energy norm. 
Rapid variations in the required neutral electronic potential in the vicinity of the nuclei are 
efficiently treated by an enriched finite element Poisson solution, using local radial solutions 
as enrichments, thus allowing a O(N) scaling solution. The formulation is presented in Sec- 
tion [2] and the solution is elaborated in Section [3j Expressions for the Coulomb energy per 
unit cell, analytically excluding the divergent nuclear self-energy, are derived. In Section HJ 
Coulomb potential and energy calculations for two canonical test cases using cubic finite 
elements (FEs) and enriched finite elements (EFEs) are presented. We demonstrate the ac- 
curacy and convergence of the approach by direct comparison to standard Ewald sums for a 
lattice of point charges, and demonstrate the accuracy in quantum mechanical calculations 
with an application to crystalline diamond. In Section [5J the main findings are summarized 
and the outlook of the proposed formulation in density-functional calculations is indicated. 



2 Formulation 

The total Coulomb potential V(x) diverges at the nuclear positions due to the divergence 
of the nuclear potentials VJ(x) = | X ^' T . of charges qi at positions Tj in the unit cell (where 
atomic units are used here and throughout). Moreover, the potentials V + (x) due to all 
nuclei and V~(x) due to all electrons in the crystal diverge at all points in the cell due 
to the long-range 1/r nature of the Coulomb interaction. To resolve these divergences, we 
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Figure 1: Schematic charge density in unit cell. Total density p = p + + p~, the sum of 
nuclear point charge and continuous electronic densities. Smooth neutralizing density p is 
introduced to render p + short-ranged and p~ neutral and amenable to direct Poisson solution, 
p is constructed as a sum of smooth, localized densities p;, strictly local within r = r c . 



introduce a smooth neutralizing density p and write the total charge p in the unit cell as 

p(x) = p + (x) + p-(x) (la) 
= p + (x) - p(x) + p-(x) + p(x) (lb) 

S v ' 11 v ' 

neutralized nuclear neutralized electronic 
charge density charge density 

= p + (x)+p-(x), (lc) 

where p + (x) = X]iPi( x ) = zCi ( ?i^( x ~~ T d is the total nuclear charge density, the sum 
of nuclear densities in the unit cell, p~(x) is the electronic charge density, and p + (x) = 
p + (x) — p(x) and p~(x) = p~(x) + p(x) are the neutralized nuclear and electronic densities, 
respectively. The integral over the unit cell J n p — = Q> the total nuclear charge in 

the cell (Fig. [I]). In order to produce a linear scaling formulation, we form p in the unit 
cell as a sum of smooth, strictly local densities pi centered at atomic positions tj with 
integrals f pi = qr- p = ^2jPi, where the sum is over all sites / in the crystal such that 
pi ^ (i.e. pi is nonvanishing) in the unit cell. Since the pi are strictly localized within a 
given radius r = r c of each site J, i.e., p/(x) = for |x — r/| > r c , the number of terms in 
the sum varies linearly with the number of atoms the unit cell. 
The total potential V(x) in the unit cell may now be written as 

V(x) = V+(x) + \/-(x) = y + (x) + V'-(x), (2) 

the sum of potentials corresponding to neutralized nuclear and electronic densities p + (x) 
and p~(x). Let V/(x) = ^1^^ and V/(x) be the potentials corresponding to charge densities 
p/(x) = g/<5(x — ti) and pr(x) at each site I and let p/(x) be spherically symmetric. Then 
both potential Vi — Vi and corresponding neutralized density pi — pi vanish beyond r = r c 
about site /, since Vi = qi/r for r > r c due to the compact-support of the associated 
spherically symmetric p/(x). The total potential in the unit cell associated with neutralized 
nuclear density p + may then be computed as a short-ranged sum in real space: 

V+(x) = J2 VKx) - V7(x) = Yl ~ KM, ( 3 ) 
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where the sum is over all sites / in the crystal such that V/(x) — V/(x) ^ in the unit cell. 
Due to the strict locality of V/(x) — V/( x ) around site /, the number of nonzero terms in the 
sum varies linearly with number of atoms in the unit cell. 

The potential associated with neutralized electronic density p~ can be obtained from a 
solution of Poisson's equation: 

V 2 F(x) = -4tt^(x), (4) 

subject to periodic boundary conditions, with continuous neutralized electronic density p~(x) 
as source term. At this point, we note that by virtue of the decomposition ([1]), the source 
term p~ in (J3J), unlike the total charge density p, is in if -1 (indeed, it is C°); so that the 
corresponding solution V~ is in H 1 , thus allowing convergent and optimal approximation 



in the energy norm ( Strang and Fix 1973 ). Moreover, the solution can be accomplished in 
O(N) operations in real space ( IBeck 20 00). where N is the number of atoms in the unit cell. 
The total all-electron Coulomb potential V(x) can thus be obtained in O(N) operations 
using decompositions <^ and (j2J), without distributed nucleus or other approximations. 

The total Coulomb energy per unit cell in the all-electron case is divergent due to the 
divergence of the nuclear self energies. Thus, the desired total Coulomb energy excluding nu- 



clear self-energy cannot be computed as in the pseudopotential case (IPask and Sterne 2005b|) 
by first computing the total Coulomb energy and then subtracting the self-energy. Instead, 
the divergent nuclear self-energy must be excluded analytically. 

The total Coulomb energy per unit cell can be expressed in terms of densities and asso- 
ciated potentials as 

E = \ I d 3 xp(x)y(x) (5a) 
Jn 

= \ [ rf 3 x(p + (x)+p"(x))(y + (x) + ^(x)) (5b) 
Jn 

= \ [ d 3 x (p + (x)V + (x) + 2p + (x)V~(x) + p-(x)V-(x)) (5c) 
Jn 

= E ++ + E + ~ + E , (5d) 

where E ++ , E + ~ , and E correspond to the p + V + , p + V~ , and p~V~ integrals, respectively, 
and the p + V~ interaction term has been retained in favor of the equivalent p~V + term to 
facilitate subsequent integration. The divergent self-energy is contained in the E ++ term. 
This may be extracted as follows: 

E ++ = \j ^p + (x)V+(x) = Ud 3 xJ2 (p/(x) - pi(x)) PA*) ~ V>(x)), (6) 
Jn Jn j j 

where the sums are over all atomic positions in the crystal with localized densities and 
potentials pi and Vj — Vj nonzero in the unit cell. Now, we restrict the neutralizing densities 
Pi to be nonoverlapping, i.e., p/(x)pj(x) = for / ^ J . Then the double summation (jBJ) 
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reduces to 



E++ = \ I E M x ) " M*)) (V^(x) - V}(x)), 

= \Y,f d ' x M x ) - *c*)) N x ) - ^w)' ( 7 ) 

the sum of neutralized nuclear self energies, where the last summation is over atoms in the 
unit cell and the integral is over all space. Extracting the divergent self-energy from ([7]), we 
have 

E ++ = E self - [ rf 3 s[p i (x)y i (x) + p i (x)(y i (x)-1/ i (x))] ) 

= E self -\Y,hVi{Ti) + / <Pxpi(x)(Vi(x) - ^(x))]. (8) 

i J 

This can be simplified further by employing a common spherically symmetric neutraliz- 
ing charge pi and corresponding potential Vi at each site. Let pi{x) = (Zi^Gx — Tj|) and 
Vi(x) = qiv(\x — Ti\), where g(r) is a smooth radial function strictly localized within r = r c , 
with J #0x1) = 1 and r c such that jOj(x) are nonoverlapping, and w(r) is the corresponding 
potential, as in Fig. [2j For r > r c then, g{r) = and v(r) = 1/r. Equation (jSJ) then becomes 

£ ++ = £ sd/ - | V [gf«(0) + q\ f ° drAnr 2 g(r)(l/r - v(r))} (9a) 
i Jo 

= E self -\Y,{q-m + qll 9 ), (9b) 

i 

where I g is the constant defined by (I9a|) . which depends only on the choice of the radial 
function g(r). In the present work, we employ a second- derivative continuous (g e C 2 (M + )) 
piecewise polynomial for g(r): 

-21(r - r c ) 3 (6r 2 + 3rr c + r2)/(57rr^), < r < r c , 
0, r > rv, 



for which 



. (9r 7 - 30r 6 r c + 28r 5 r 2 - Ur 2 r 5 c + 12rJ) /(5rj!), < r < r c , 
1/r, r > r c 



c- 



and 

7 5 = 10976/(17875 r c ). (10c) 
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Figure 2: Smooth, strictly local unit charge g(r), corresponding potential v(r), and 1/r, with 
cutoff radius r c = f . g(r) = and v(r) = 1/r for r > r c . 

From (j5cl) and (]5d[) . the interaction term is given by 

E + - = [ d 3 xp + (x)V~(x) 
Jn 

= [ rf 3 x^( P/ (x)-p/(x))^(x) 
Jn j 

= J2^( T ^- I d s xp(x)V-(x), (11) 
i J n 

where the sum on / extends over all atoms in the crystal with p~j ^ in the unit cell and 
the sum on % is over atoms in the unit cell. 

Being free of nuclear singularities, the E term can be evaluated straightforwardly as 

E— = \ [ Ap"(x)r(x). (12) 
Jn 

Collecting E ++ , E^ , and E from (19bj) . (TTTT) . and (fT2j) above, we arrive at the following 
expression for the all-electron Coulomb energy per unit cell, excluding nuclear self-energy: 

E-E self = [^-(r^-^HO) + /,)]+ f rf 3 x(ip-(x)-p(x))V-(x), (13) 
. Jn 

where the integral is over the unit cell and the sum is over atomic positions in the cell. Since 
the densities and potentials in (fT3"j) can be obtained in O(N) operations, as described above, 
the energy too, as formulated in (1131) . can be obtained in O(N) operations. 

In the above expression, the p + V~ interaction term was retained in favor of the equiva- 
lent p~V + term to facilitate analytic integration. However, whereas this eliminates a three- 
dimensional numerical integration, it produces a term in the energy requiring a pointwise 
evaluation of the Poisson solution V~ . In basis-oriented, variational approaches for the Pois- 
son solution, however, energy integrals of the solution converge more rapidly than pointwise 
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values QStrang and Fix 1973[ ), and so it is of interest to develop an alternative expression 
free of pointwise evaluations of V~ . This can be accomplished at the cost of an additional 
numerical integration by formulating the E A interaction term based on the corresponding 
p~V + integral rather than p + V~ integral as in (l5cj) . 
The interaction term is 

E + - = [ d 3 xp + (x)V-(x) = [ d 3 xp~(x)y+(x). (14) 
Jn Jn 

This may be expressed as 

e + -= [ dMp-(x) + £Mx)]E(^( x )-^( x )) 

Jn i j 

= W rf 3 xp-(x)(V 7 (x)-V}(x))+E / rf 3 xp 7 (x)(V 7 (x)-V7(x)) 
j Jn j Jn 

= £/ rf 3 xp^(x)(V^(x) - V^x)) + [ <Pxpi(x)(yi(x) - K(x)) (15) 

i i 

for nonoverlapping densities p/(x), where / extends over all contributing sites in the crystal; 
i, over positions in the unit cell, and the final integrals are over all space. On letting 
p~i(x) = qig{\x. — Ti\) and Vj(x) = gjw(|x — Tj|), where g(r) and v (r) are as in ( 110 aft and 
(llObp . the interaction term becomes 



E + - = Y,^ [ d 3 xp'^)(l/r l -v(r l ))+J2'l"l g ^ 



(16) 



where is the sphere of radius r = r c within which density p~j(x) is localized, = |x — Tj|, 
and I g is as in (llOcp . With this formulation of the interaction term, the all-electron Coulomb 
energy per unit cell, excluding nuclear self energy becomes 

E - E S(df = J2±q?(l g - v{0)) +J2q t [ ' d 3 x p-(x)(l/r, - v(n)) +± [d 3 x pT {x)V~ (x). (17) 
i i J^ J Q 

In the special case of constant p~(x), as in the Ewald problem, the integral in the second 
term of (fT7|) reduces to I sp h = 14p~7rr^/75. In the more general case, the weak angular de- 
pendence of p~ (x) about each site % allows for efficient evaluation using Gaussian quadrature 



in spherical coordinates (IStroud 1971[) . 



3 Solution 

The computation of the all-electron Coulomb potential and energy as formulated in (J2J), (jl3j) . 
and ( 1T7|) requires the solution of the Poisson equation (J3J) for the potential V~ corresponding 
to the neutralized electronic charge density p~, which may be accomplished by a number of 



methods. Here, we employ an enriched finite element (EFE) method (Strang and Fix 1973 



IMelenk and Babuska 1996; Babus ka and Melenk 1997] ISukumar and Pask 2009)) in order to 



efficiently resolve the sharp variations in the all-electron densities and potentials while re- 
taining systematic improvability and O(N) scaling in the solution process. Note that some 
form of multilevel preconditioning is generally required to achieve linear scaling of solution 
time with respect to number of degrees of freedom (DOFs) in FE and related methods. 
While not demonstrated as yet in the context of EFE methods, to our knowledge, it is ex- 
pected that such preconditioning will apply here also. EFE methods add fixed enrichment 
functions to the classical FE basis in order to efficiently represent rapid variations in the 
solution. This leaves a smooth residual (difference between exact solution and enrichment) 
for the remaining FE basis to represent, thus allowing a substantially coarser FE mesh. EFE 
methods can thus be understood as FE methods on softer problems, and so obtain classical 
asymptotic convergence rates with respect to mesh size, depending only on the completeness 
of the classical part of the basis (as shown below). Hence, multilevel preconditioning may 
be expected to apply to EFE as for FE. Indeed, consistent with this expectation, linear 
scaling with respect to number of DOFs has been demonstrated recently in the context of 
partition-of- unity enrichment methods ( ISchweitzer 20031 ISchweitzer 2008|) . 

In the classical FE Poisson solution (IPask and Sterne 2005a|) . the potential V~ is ex- 
pressed as a linear combination of strictly local, piecewise polynomial basis functions {<pi}: 

Vf E (x) = J2^ a ? E - ( 18 ) 

i 

In the enriched FE solution, a set of functions {ip a } which incorporates a priori information 
about the solution is added to the classical basis in order to substantially reduce the number 
of basis functions required to attain a given accuracy; so that V~ is expressed as 

y-(x) = ^0 l (x)a i + ^^ a (x)6 a = ^$ fc (x)c fc , (19) 

i a k 

where a is summed over the atoms and = {<pi} U {ipa} is the enriched FE basis, the 

combined set of classical FE and enrichment basis functions. 

In the present case, the desired solution V~ = V~ + V is the potential associated with 
neutralized electronic charge density p~ = p~ + p, where p = YliP~i> the sum of smooth, 
strictly local, nonoverlapping neutralizing densities over atomic sites I in the crystal such 
that pi ^ in the unit cell. In the all-electron case, the electronic density p~ is large 
in magnitude, rapidly varying, and isolated-atom-like in the vicinity of the nuclei while 
much smaller in magnitude and more smoothly varying elsewhere. The neutralizing density 
p is moderate in magnitude and smoothly varying throughout the cell. Hence, the rapid 
variations in the desired solution V~ associated with neutralized electronic density p~ are 
confined to the vicinity of the nuclei where the electronic density p~ is large, rapidly varying, 
and atomic-like. To increase the efficiency of the representation, therefore, we might add to 
the basis the potentials vj corresponding to the neutralized electronic densities pj = pj + pi 
in the vicinity of each site /; where the local electronic densities pj vary like the total 
electronic density p~ in the vicinity of site / and integrate to the appropriate charge — qj. 
These may be obtained, for example, from a partitioning of the self-consistent electronic 
density or from isolated atom densities. 
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For each atom a in the unit cell, then, an enrichment function 



\C(x) = 5Jf;-(x-R) (20) 

R 

approximating the desired solution V~ in the vicinity of atom a is constructed, where R 
denotes lattice translation vectors such that 0~(x — R) ^ in the unit cell. Since the 
enrichment is only needed in the vicinity of atom a, the v~ and associated summation 
are short-ranged. The enriched basis functions ip a in ( fl9|) are taken, then, as the periodic 
enrichment functions V~: 

^(x) = V7(x). (21) 

To solve the Poisson equation subject to periodic boundary conditions, it is sufficient 
that the basis {$£,•} satisfy 

$ fe (x + R) = $ fc (x) (22) 
for all x on the boundary (IPask et al. 19991 IPask and Sterne 2005a|) . Classical basis func- 



tions {<pi} satisfying this condition are readily constructed (IPask et al. 1999]) . Furthermore, 
the enriched basis functions {^ a } in ( )2~Tj) are periodic in the unit cell by construction 
(Eq. (12TJI)). Thus the enriched FE basis {&k} as a whole satisfies the required condition 
also ( ISukumar and Pask 20090 . 

Having so constructed the enriched basis satisfying the required boundary con- 

ditions, the enriched FE solution V^ _ (x) = Ylj ^j( x ) c j °f the Poisson equation (j3J) with 
neutralized electronic density p~ (x) is then determined by the sparse, symmetric linear sys- 
tem (IPask and Sterne 2005aP 



jCj — fi, (23a) 



where 



Lij = / AV$i(x) • V$j(x), (23b) 
Jn 

fi = 4n f c/ 3 x$,(x)p-(x). (23c) 
Jn 



4 Results 



We demonstrate the accuracy and convergence of the formulation presented in Section [3] 
on two canonical test cases: the Ewald energy of a bcc crystal, for which a reference value 
is available via standard Ewald summation (lEwald 19211) . and the all-electron Coulomb 



potential and energy of crystalline diamond. In the numerical computations, the paral- 
lelepiped unit cell is discretized into regular m x m x m serendipity cubic (32-node) finite 
elements. The shape function expressions for serendipity cubic brick elements are provided 
in Zienkiewicz et al. (2005) Upon obtaining radial atomic potential solutions v~ analyt- 
ically or from a one- dimensional solver (exploiting spherical symmetry), the neutralized 
atomic potential functions v~ are tabulated on a one-dimensional equi-spaced grid for the 



11 



Ewald problem and on a logarithmic grid for the diamond problem. A quintic spline-fit of 
the tabulated data is formed and the enrichment function is constructed as in (1201) . which is 
then used in the numerical computations. Numerical integration for the FE computations 
is carried out using 5x5x5 Gauss quadrature, and higher-order quadrature is employed in 
EFE computations. The FE computation on a m x m x m mesh has 7m 3 degrees of freedom, 
and the corresponding EFE computation has just two (number of atoms) more. 

4.1 Ewald problem 

For the classical Ewald energy, we consider a reference bcc crystal with unit atomic spacing 
defined by lattice vectors 

ai = a(l,0,0), 
a 2 = a(0,l,0), 

33 = 0(0,0,1), 

with unit point charges g« = 1 located at positions (in lattice coordinates) 

n = (o,o,o), 

r 2 = (1/2, 1/2, 1/2) 

and a constant negative background p~ = — 2/a 3 bohr -3 , where a = 2/y3 bohr. The Ewald 
energy computed by standard N 2 scaling summation is —1.5758343085 Hartree/atom. 

To compute the energy using the O(N) formulation (fT3"l) or (fT7|) . we take g(r) as in (j 10 aft 
with r c = 1/2 bohr, so that the resulting neutralizing densities p~/(x) = p(|x — tj|) at each 
site I in the crystal are as smooth as possible without overlapping. The total neutralizing 
density in the cell is then p(x) = ^jPj(x), where the sum over positions I in the unit and 
nearest neighboring cells is sufficient due to the short-range of g(r). For the sake of clarity, 
we distinguish between the cut-off radius r c used to represent the neutralizing charge pi and 
that used to form the enrichment function corresponding to the neutralized electronic charge 
density. We refer to the former as r cn and to the latter as r ce . Since the electronic density 
p~ is constant, we take as local electronic densities p7( x ) = — 9(\ x ~ T A) with r ce = 1 a.u. 
so that Yli Pj approximates the constant p~ in the cell with I running again over positions 
in the unit cell as well as the nearest neighboring cells. The potentials corresponding to 
neutralized electronic densities pj = pj + pi at each site I are then vj = vj + vi, where 
vj(x) = — v (|x — Tj\) with r ce = 1 a.u., £j(x) = u(|x — Tj\) with r cn = 1/2 a.u., and v(r) is 
as in (llQbl) . The enrichment functions V~ are then formed as in (|20|) . The required potential 
V~ corresponding to neutralized electronic density p~ = p~ +p is obtained from the enriched 
FE solution of the associated Poisson equation subject to periodic boundary conditions, with 
enrichment functions V~ . The total (electronic + nuclear) potential is then obtained from 
(j2J) and the energy, from ffTB"]) or ffTTj) . 

The numerical results for the Ewald problem are shown in Figure [31 The electronic, 
neutralizing, and neutralized electronic charge densities are shown in Figure Eh- The FE 
potential solution corresponding to the neutralized electronic charge density is plotted along 
the diagonal in Figure [3b for different meshes; convergence is observed as the number of 
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elements along each coordinate direction increases from 3 to 5 to 7. The two enrichment 
functions, one for each atom, are shown in Figure |3b, and the EFE solution for a 12 x 12 x 12 
mesh is depicted in Figure |3Ji- The total potential solution, which includes the singular nu- 
clear potential contributions, is shown in Figure |3fe. Figure shows the convergence curves 
for the Coulomb energy per atom for FE and EFE solutions. Numerical integration in the 
EFE solution is carried out with a 20 x 20 x 20 Gauss quadrature rule in each element, which 
ensures that the quadrature error is less than the approximation error. From Figure H^, we 
see that use of the pointwise expression f[T3~j) adversely affects both accuracy and rate of con- 
vergence. To explore this further, the pointwise error Yfi=i V~( T i) * s plotted in Figure Hb, 
where the EFE solution (16 x 16 x 16 mesh) is used as reference. Clearly, the pointwise error 
is appreciable, which explains the decrease in accuracy and the non-monotonic convergence 
in the curves plotted in Figure UK The enriched finite element solution provides an accuracy 
of order 10~ 8 Ha in the Coulomb energy with 9 elements in each direction whereas with finite 
elements, 32 elements in each direction is required for the same accuracy. On using the inte- 
gral expression in ( fTT|) to compute the Coulomb energy, the optimal sextic rate of convergence 



is obtained with FE and EFE (Figure H^i), consistent with theory (Strang and Fix 1973). 

From the error curve for EFE in Figure (integral expression) , we note that the mesh 
4x4x4 delivers much better accuracy and convergence rate than the other meshes. Further 
analysis helps to explain this anomaly. For the Ewald problem, the corner atom is at the 
origin and the center atom is at r = (1/2,1/2,1/2) (length of the diagonal is 2 a.u.), 
and r cn = 1/2 a.u. ensures no overlap between the neutralizing densities of the two atoms. 
Furthermore, the support of the neutralizing charges coincides with the location of the nodes 
(4x4x4 mesh) along the diagonal of the cube. We repeat the calculations for two more 
cases: r cn = 1/3 a.u., r ce = 1 a.u.; and r cn = 1/3 a.u., r ce = 2/3 a.u. The results obtained 
for all three cases are shown in Figure |5j One can observe that the solutions for the mesh 
4x4x4 with r cn = 1/3 a.u. are consistent with the solutions on other meshes, and hence 
the markedly better accuracy for the case when r cn = 1/2 a.u. (Figures H^i and [5]) is a 
consequence of the choice r cn = 1/2 a.u. 

Numerical integration to compute the weak form integrals is carried out using standard 
tensor-product Gauss quadrature. Since the right-hand-side of the Poisson equation (jlj) is 
not a polynomial, the degree of the quadrature rule must be selected so that the integration 
error is at least an order smaller than the approximation error in order to obtain variational 
results. Figure E] shows the error in the Coulomb energy per atom (24 x 24 x 24 FE mesh) as 
a function of the number of quadrature points in each direction. The error in the FE solution 
is shown by the horizontal lines, and Figure |6] reveals that at least five points (5x5x5 
quadrature rule) are needed to ensure that quadrature error is below approximation error. 

4.2 Diamond 

We now consider diamond, with unit cell defined by lattice vectors 

ai = a/2(0,l,l), 
a 2 = a/2(l,0,l), 

33 = 0/2(1,1,0) 
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Distance along diagonal (a.u.) Distance along diagonal (a.u.) 



(c) (d) 




Figure 3: Finite element (FE) and enriched finite element (EFE) solutions for the Ewald 
problem, (a) Charge densities; (b) FE solutions; (c) Enrichment functions; (d) EFE solution 
(12 x 12 x 12 mesh); and (e) Total EFE solution (12 x 12 x 12 mesh). 
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Figure 4: Convergence curves, (a) Error in Coulomb energy per atom. Convergence rates: 
FE (integral) is 6.06, EFE (integral) is 6.15, FE (pointwise) is 4.97, and EFE (pointwise) is 
3.83. (b) Error in Y^=i V~{ T i)- EFE solution on 16 x 16 x 16 mesh is used as reference. 
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Figure 5: Error in Coulomb energy per atom for EFE computations using different r cn and 



and carbon atoms at positions 

n = (0,0,0), 

r 2 = (1/4, 1/4, 1/4), 

where a = 6.75 bohr. For the all-electron problem, we consider a total charge density p 
consisting of nuclear point charges qt = 6 at positions Tj and electronic density p~ from 
overlapping all-electron atomic densities. 

To compute the energy using the formulation ( {TBI or i fTTl) . we take g(r) as in f llOap with 
r c = 1.4 a.u., so that the neutralizing densities p~/(x) = 6g(|x — tj|) are nonoverlapping. 
The total neutralizing density in the cell is p(x) = X]/Ar( x )' where the sum over positions 
I in the unit and nearest neighboring cells is sufficient. The electronic charge density is 
p~ = ^2 1 Pj , which is the sum of all-electron atomic densities, where the sum over positions 
in nearest and second-nearest neighbor cells is sufficient to capture the tails of the pj. The 
potentials corresponding to neutralized electronic densities pj = pj + pj at each site / are 
then vj = vj + Vi, where vj is the all-electron atomic potential and vi is the potential 
corresponding to pi. The integral in the second term of ( 1T7|) involving p _ (x) is computed 
within a sphere of radius r c using Gauss quadrature in spherical coordinates. A tensor- 
product quadrature rule with 531 points (spline-fit has 177 knots) in the radial direction, 
and 20 points in each of the two angular directions is used. We obtain I sp h = —9.44150186 
with precision to all digits shown. The enrichment functions V~ are formed as in (120]) . 
with the vj brought smoothly to zero at radius a/y/3 in order to maximize their extents 
consistent with summation over unit and nearest neighbor cells only. The required potential 
V~ corresponding to the neutralized electronic density p~ is then obtained from the enriched 
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Figure 6: Comparison of quadrature errors for FE computations on 24 x 24 x 24 mesh. 
Horizontal lines show the best solution error obtained for each case. 



FE solution of the associated Poisson equation, with enrichment functions V~; whereupon 
the Coulomb potential and energy are obtained from (J2J) and ([TBI) or (fT7|) . respectively. 

The electronic, neutralizing, and neutralized electronic charge densities are shown in Fig- 
ure Uh- The neutralized finite element potential solutions along the diagonal of the unit cell 
are plotted in Figure [7b for 8, 16, and 32 elements along each primitive lattice direction. The 
two enrichment functions, one for each atom, are shown in Figure [7b, and the neutralized 
EFE potential solution and total potential solution for the 24 x 24 x 24 mesh are illustrated 
in Figures [TH and [Tfe. The error in the Coulomb energy with mesh refinement is plotted 
in Figure [7f for FE and EFE methods; with EFE result on a 24 x 24 x 24 mesh as reference. 
The enriched finite element solution has an accuracy of 4 x 10 -6 Ha in the Coulomb energy 
for a 16 x 16 x 16 mesh (28,674 degrees of freedom), whereas the the best uniform-mesh FE 
result provides an accuracy of only 6 x 10 -3 Ha on a 32 x 32 x 32 mesh (229,376 degrees 
of freedom). While the use of adaptive higher-order finite elements will require fewer basis 
functions than uniform FE, previous studies (IGoedecker and Ivanov 19981 ITorsti et al. 2006"! 
Bylaska et al. 2009} ILehtovaara et al. 2009[ |Suryanarayana et al. 2010[) suggest that their 



performance in terms of number of basis functions for a prescribed accuracy may not com- 
pare favorably to EFE since they do not incorporate physics-based knowledge (nature of 
variations, in addition to scale) within the approximation space. 

The numerical integration for FE and EFE solutions is performed using tensor-product 
Gaussian quadrature rules. For EFE computations, the number of integration points are 
increased until the integration error is below the approximation error. Due to the large 
values and sharp variations of the enrichment functions in the vicinity of the atomic posi- 
tions, higher-order quadrature is required in finite elements that are close to the atoms. As 
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Figure 7: Finite element (FE) and enriched finite element (EFE) solutions for crystalline 
diamond, (a) Charge densities; (b) FE solutions; (c) Enrichment functions; (d) EFE solution 
(24 x 24 x 24 mesh); (e) Total EFE solution (24 x 24 x 24 mesh); and (f) Error in Coulomb 
energy per atom (integral expression). Convergence rate for EFE: meshes 8-12-16 give 5.38 
and meshes 12 and 16 yield 5.79. EFE solution on 24 x 24 x 24 mesh is used as reference. 
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a result, we use higher-order quadrature in elements with at least one vertex at an atomic 
position, and a lower-order quadrature rule in all other elements. The reference solution 
for the Coulomb energy per atom in Figure [7f is —75.985203 Ha, which is obtained using 
a 25 x 25 x 25 quadrature rule in elements that are in the vicinity of the atoms and a 
15 x 15 x 15 quadrature rule in other elements. The integration scheme we adopt is the 
simplest rule that provides us with the required accuracy. For better efficiency, a possi- 
ble approach is to use tetrahedral mesh generation techniques to construct a partition of 
a finite element with an atom (Coulomb singularity) located at one of its vertices or in its 
interior. With a graded tetrahedral mesh that is focused towards the atomic position, stan- 
dard quadrature rules within the tetrahedral elements would suffice to accurately integrate 
the neutralized charge density and the enrichment functions. Within an EFE method, the 
tetrahedral mesh so constructed will be solely used for the purpose of numerical integra- 
tion; the number of degrees of freedom remain unchanged. For the Coulomb singularity, 



Batcho (2000) and Havu et al. (2004) adopt the Duffy transformation (Duffy 1982) to nu- 



merically integrate matrix-elements with 1/r factors. The development of accurate and 
efficient quadrature schemes in EFE methods is a topical issue at the forefront of current re- 
search; recent work in this direction includes a generalization of the Duffy transformation for 
integrating power singularities ( IMousavi and Sukumar 2010|) . Adaptive integration schemes 
or quadrature rules that are tailored to the form of the enriched FE basis functions could 
prove to be more efficient to solve the neutralized Poisson problem. This and related topics 
are subjects of ongoing research. 



5 Conclusions 

We have presented a linear scaling formulation for the solution of the all-electron Coulomb 
problem in crystalline solids. The resulting method includes full nuclear potentials, with 
no smearing approximations, is systematically improvable, and well suited to large-scale 
quantum mechanical calculations in condensed matter in which the Coulomb potential and 
energy of a continuous electronic and singular nuclear density are required. Linear scal- 
ing is achieved by introducing smooth, strictly local neutralizing densities to render nuclear 
interactions strictly local, and solving the remaining neutral Poisson problem for the elec- 
trons in real space. Rapid variations of the electronic density in the vicinity of the nuclei 
were efficiently treated using enriched FE methods, with isolated atomic solutions as enrich- 
ments. By considering different interaction terms, we derived two equivalent expressions for 
the Coulomb energy per unit cell — one involving a pointwise evaluation of the neutralized 
electronic potential, the other requiring the evaluation of a spherical integral. By avoid- 
ing pointwise evaluation of the C° FE solution, the integral expression proved superior in 
both accuracy and convergence rate with respect to the number of elements. For the Ewald 
problem, accuracy of order 10~ 8 Ha was realized with enriched FE on a mesh with 5105 
degrees of freedom. Comparable accuracy on a uniform FE mesh required 229,376 degrees of 
freedom. For the Coulomb energy of diamond, the enriched FE solution with 28,674 degrees 
of freedom attained an accuracy of 4 x 10~ 6 Ha. 

While the calculations here employed a finite element basis, the formulation applies quite 
generally to any desired basis for the residual Poisson solution and so should prove of wide 
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applicability. For example, in the context of high-accuracy Gaussian based calculations, 
where fast multipole methods become more costly, the present method may provide an at- 
tractive alternative for O(N) parallel solution. In such a case, the strictly local polynomial 
neutralizing functions employed here might be replaced by correspondingly localized Gaus- 
sians. Any convenient functional form satisfying the required non-overlap conditions can be 
employed to equal effect. 
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